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IMPLICIT APPROXIMATE- FACTORIZATION SCHEMES FOR THE 
LOW-FREQUENCY TRANSONIC EQUATION 
William F. B^^llhaus and Joseph L. Steger 

Ames Research Center 
and 

U.S. Army Air Mobility R&D Laboratory 
Moffett Field, California 94035 

SUMMARY 

Two- and three-level implicit finite-difference algorithms for the low- 
frequency transonic small disturbance-equation are constructed using approxi- 
mate factorization techniques. The schemes are unconditionally stable for the 
model linear problem. For nonlinear mixed flows, the schemes maintain stabil- 
ity by the use of conservatively switched difference operators for which sta- 
bility is maintained only if shock propagation is restricted to be less than 
one spatial grid p' it per time step. The shock-capturing properties of the 
schemes are studii r various shock motions that might be encountered in 
problems of engineering interest. Computed results for a model airfoil prob- 
lem that produces a flow field similar to that about a helicopter rotor in 
forward flight show the development of a shock wave and its subsequent propa- 
gation upstream off the front of the airfoil. 


I . INTRODUCTION 


Frequently, one employs simplified equations that extract the essential 
physics from a more corplete set of equations that model some physical phenom- 
enon. These simplifications, while limiting the range of applicability, often 
significantly reduce the complexity and expense of computing solutionc. Wlien 
such approximations are made, the numerical algorithm used to solve the sim- 
plified equations should be matched to the same range of applicability. 

Of interest here is the solution of the transonic flow field about a thin 
airfoil executing some arbitrary small amplitude oscillatory motion in a uni- 
form stream. The small disturbance and transonic assumptions allow simplifi- 
cation of the governing gas-dynamic equations to (ref. 1, p. 4) 

k2M„2<J.,t + 2kM„24.^^ = + ,|,yy (1) 

where = 1 - _ (y + ij) is the disturbance velocity potential, 

Y is the ratio of specific heats, and m is a function of (see refs. 2 
and 3). The coordinate system is fixed with respect to the airfoil, and x 
is alined with the free-stream direction. The flow is defined to be locally 
subsonic or supersonic, relative to the fixed coordinate system, for V,, > 0 
or V(. < 0, respectively. A measure of the degree of unsteadiness of the 
motion is given by the reduced frequency k; for an airfoil of chord length 
' , with free-stream velocity U„, and executing some periodic motion of 
frequency m, 
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( 2 ) 


(i)C 


The variables <t>, t, y, and x have been nondlmenslonallzed by cU„, 1/ai, c, 
and c, respectively. In the derivation of (1), it is assumed that the flow 
is inviscid and transonic (M„(, 1) and that velocity disturbances due to the 

airfoil are small compared to U„. Further simplification can be achieved by 
restricting the range of reduced frequencies, as we shall see. 

The traces of the characteristic surfaces in tl'S plane of the airfoil for 
the linear form of (1) are shown in figure 1. Waves that advance in the flow 
direction propagate downstream relative to the airfoil with velocity 
(Mcb + l)/Mt„k in the nondiraensiuiial coordinate system. Waves that recede 
(relative to particle paths) travel upstream against the flow with velocity 
PluLd particles travel with velocity 1/k (i.e. , Uo, in the 
dimensional system). For low-frequency transonic flows (ref. A), 

k ~ 1 - ~ << 1 (3) 

where t is the airfoil thickness-to-chord ratio. So dist bances traveling 
along the advancing waves propagate rapidly away from the air roil, while those 
traveling along the receding waves remain closer to the airfoil and dominate 
the solution there. For low-frequency motion, (1) reduces to the low-frequency, 
transonic, small-disturbance equation 

2kM„'4>xt = + -i-yy 

which is subject to the restrictions (3). The numerical solution of (A) is 
useful, for example, in aerodynamic flutter analysis and the treatment of high 
speed helicopter rotors. 

An efficient computational algorithm for either (1) or (A) should have a 
time step restriction for stability that is no more restrictive than the time 
step required to adequately resolve the unsteady flow field. Since, for low- 
frequency motions, the solution is dominated by disturbances traveling along 
the receding waves, the time step should be related to the time scale of 
receding wave propagation, as well as the wave length of these disturbances. 
Explicit schemes for transonic flows have time step restrictions for stability 
that are based on the advancing wave time scale. Hence, for low-frequency 
motion, they have time step restrictions for stability that can be orders of 
magnitude more severe than required for accuracy. This problem of stiffness 
was overcome in reference A by the use of semi-implicit difference schemes 
applied to (1) and (A), with a resulting time step restriction based solely cn 
the receding wave time scale. 

Results obtained using the semi- implicit scheme applied to (A) have been 
compared with solutions to the complete Eulerian gas-dynamic equations in 
references 2 and 3. The solutions compared favorably for cases in which the 
assumptions (3) were valid, and the simplified theory resulted in a consider- 
able savings in computer run time and storage. However, it will be shown in 
what follows that considerably larger time steps than allowed by the semi- 
implicit scheme stability requirement 
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2kM„2^5j 

At i 22 ^ 

1 - M„2 - (Y + 1)M«,%| 


( 5 ) 


can be taken without sacrificing accuracy for truly low-frequency notions. 

Note that (5) is especially restrictive near the airfoil leading and trailing 
edge singularities, where |4>^| is large. 

Two fully implicit schemes are reported here that have no time step 
restriction based on a linear stability analysis. One scheme, an alternating- 
direction implicit algorithtif, is applied to (4) in a manner similar to the 
Douglas-Gunn algorithm for the heat equation (ref. 5). The other scheme, an 
apprrximate faccorizatlon scheme, also factors difference approximations to 
(4) such that simple matrix solutions are required first in one spatial direc- 
tion and then the other. In Section II, both schemes are investigated in 
terms of stability and accuracy for a model linear equation. In Section III, 
the shock capturing characteristics of the schemes are investigated, and con- 
ditions to preve' t nonlinear instability and shock-capturing overshoots are 
found. In Section IV, the difference schemes are presented for the complete 
equation (4). Finally, in Section V the algorithms are applied and results 
compared for the simulated motion of a helicopter blade tip. 


II. ALGORITHMS FOR THE MODEL LINEAR EQUATIONS 


We begin by considering the model linear equation for (4) 

*^xt ° ^'^xx ‘f’yy 

An implicit differencing of (6) is sought to remove any time step restrictions 
other than those imposed by considerations of accuracy, and, for reasons of 
efficiency, computer-time consuming inversions must be avoided. The system of 
equations generated by implicit difference operators is perhaps most effi- 
ciently solved by structuring the difference operators so that they can be 
factored into easily inverted products, For simple difference equations one 
can sometimes construct fast direct solvers in this way, but £s the complexity 
of the equations increases, approximate factorizations are usually sought. 

This is the conropt of splitting and alternating directions (see, e.g., 
refs. 5, 6, and 7) which has recently been applied to the Eulerian gasdynamic 
equations (ref. 8). (The terminology "approximate factorization" is adopted 
from Yanenko (ref. 6, p. 27) and is perhaps the most descriptive name for this 
general class of procedures.) In this section two such procedures are devel- 
oped which shall be referred to as a three-level approximate factorization 
(AF) scheme and a two-level altetmating direction (ADI) scheme. 


a) Approximate Factorization Algorithm 


An implicit differencing of (6) for 6 > 0 is 

„ ,n+l o- ,n+l . - 


(7) 
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where 


(AxAt)"\l - E'l)(l - E-1) , 

(A..)-2(e"1 - 2 + eJI) , 

(Ay)“2(E-l - 2 + E+J) . 


"xt 


"XX 


'yy '•—J ' '-^y “ “y 

Here the shifting operators are defined by 


[0(At) + 0(Ax)] 
[0(Ax^)] 

[0(A/2)] 


±1 n nil 

E 6 , « 4 

t ‘Pj.k ^j,k 


:tl n _ n 
®x *^j,k ■ '**j±l,k 


(8a) 

(8b) 

(8c) 


(9) 




= if 


n 

J.k±l 


with t = nAt, X = jAx, and y = kAy. In the above, kept first-order 

accurate for simplicity — a second-order-accurate version will be described 
later. Note that (8b) can be rewritten 6^^ = (Ax)“2(e+^ - 1)(1 - E“^). 


Various factorizations of (7) are pofsible, but one that makes use of the 
common factor in 6jjx 6^t requires only simple bidiagonal and tri- 

diagonal inversions is 

[1 - a(4‘ - 1)][1 - e;‘ - bSyj,]*"';-' . Cl - e; 1 + ab(E+'- 1) (10) 

where a = S and b = AtAx. The approximate factorization (10) differs 
from (7) by the added cross term 


[ab(E^^ - l)6yy][E^^ - = 0(AxAt3<t>xyyt) 

which, compared to the other terras, is of second order. Further, if a steady 
state is reached, the added cross term is zero. 


The advantage of the AF form is that the inversion process is quite 
simple. As indicated in sketch (a) only bidiagonal and tridiagonal inversions 
are required, and the intermediate solution results can be overloaded into the 
same double array. A difficulty with the AF method is the need to supply 
the value (1 - - b6yy)(J)";J;l as a boundary condition during the inversion 

of each upper bidiagonal. For example, if given on the boundary, 

^n+1 remains unknown, unless one can estimate | on the downstream boundary 
' J Ic * 

A simple approach is to use data at level n in place of the unknown n+1 
boundary data. However, this particular approximation can impose conditional 
stability criteria, although not very stringent ones, especially if the local 
value of Ax is large. 
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Solution Procedure for equatirn (10) using a double array, ij)” 


Form right-hand side (overstore . Into , ) 

J f K ' ^ K 


,j+i 


Inversion of first factor (overstote . into iji ) 

J » ^ J 


♦j.k’ <*j,k + “Vi-k''” * 


j J.l 


with 


.n +1 


•^J+1 ,k 




Tridiagonal inversion (overstore into ij)" , ) 

j * j » 


<!> 


n+l 


n +1 


j,k (Ay)2 ^‘*‘j,k-l " ^^j,k “^j.k+l^ " “^j.k ’*’j-l,k 


k = 1,K 


i *h 


1 ,K 


J = l.J 


Sketch (a) 


The algorithm (10) is unconditionally stable for S > 0 by a linear sta- 
bility test (see Appendix A, part a, for details). If 8 < 0, the algorithm is 
modified so that 6 ,^ 3 ^ is the first-order three point backward difference 
operator, and stability is again unconditional. The AF algorithm for 8 < 0 
is given by 

[1 - a(l - e; 1)][1 - - b6yy]4."'J'J. = [1 - e ;1 + ab(l - E;l) 6 yy](|»"^j^ (12) 

and the solution process is much as in sketch (a) — a lower-bidiagonal inver- 
sion replaces the upper-bidiagonal inversion, and boundary conditions are 
simplified. 

Second-order accuracy in time can be obtained by using the following 
three level finite-difference operator for 

63 , 4 . = (2AtAx)“l[3- Et^(2- Ex)](l- E"^)(l- E“l) , [O(At^) + 0(Ax2)] (13) 

The operator (13) was derived by Taylor series expansion so that only the fac- 
tor (1 - E^^) appears at the n +1 level, and the factorization is much as 
before. The three-lo.vel scheme can be programmed so that only two levels of 
computer storage are required. 

The extension of the AF algorithm to three dimensions is easily accom- 
plished. For the equation 
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( 14 ) 


“^yy + 'I’zz 

an unconditionally stable approximate factorization using operators (8) with 
6gjj similarly defined is given by 

[1 - a(E+l - l)][(l - b6yy)(l - b6^^) - 

= [1 - E;;! + a(E+l - l)b(6yy + - bfiyy6„) b2fiyy6^j4.j^^ (l5) 

The inversions are again straightforward, although a complicated cross term 
has to be computed at level n. In many problems the term b26yy62j could be 
safely neglected and always without losing second-order accuracy. 


b) The Alternating-Direction Algorithm 

A two time level differencing f (6) is given by 


^xt'*’ 


ii+i _ _1^ f . 

j,k " Z 


»yy> + 




(16) 


where the previous operators (8) are used; or, to maintain second-order 
accuracy 


'xt 


= (AxAt)"^(l - EI^) 


- F"l' 


h - E-‘1 




(1 - E'^) , [0(At2) + 0(/x2)] (17) 


An approximate factorization of (16) is obtained from the Douglas-Gunn 
alternating-direction algorithm with the identity operator replaced by as 

follows 


R At i 

■ 2 ^yy 


ffi.. - 


) 


[^x-TS'xxJ^J.k 


•rn+l 


, n+1 


*^x 2 ^yy ^^^'^xx)'}’j 

^x-(>j.k - T 


where 


or 


= (1 - e'S/Ax 


(first order) 


3 - e;;> 

1 - e;"' 

2 

V < 

1 


(second order) 


(18a) 

(18b) 


The first step of this procedure is a consistent approximation to (6) and, in 
fact, is identical to the semi- imp lie it method of reference 4. 
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To show chat (18) represents an approximate factorization to (16), Insert 
(the constant of summation is taken as zero) between [6v“ (6t/2)^w] 
and 4*. 1, in (18a), and then eliminate & 4 , using (18b) with the result 

J ,x X J ,K 


(6, 


X - T «yy>‘x‘[<«x - f 6‘xx)*”i + f ®‘xx'»j.kl 


At 


n+l . At 


“ (<Sx +T Sy ^‘^s^xx)<!>5,k 


Algebraically rearranging terms 



At ^ At 

2 ®yy ■ 2 


66 


xx^ 


n+l 

J.h 




j»h 



^'^yy^x^'^xx(^ 



)<t> 


n+l 


( 20 ) 


With At divided out the error term is 0(ac^<pyyxt^ ns for the AF scheme. 

n4"l 

In the inversion procedure for (18), 4 must be supplied as a boundary 

condition in step (18a). From step (18b) this means that 



n+l At ,_i. /.n+l .n . 

'j,k-T''x''xx(*j,k- +j,k> 


must be supplied on the boundaries, but to second-order accuracy 4j,k be 

taken as k °n the boundary. This is typical of the Douglas-Gunn metbud — 
implementation of other, even rearranged, approximate factorizations may 
require that the function and a derivative be specified on the boundary. 


For 6 > 0, (18) is unconditionally stable by a linear stability test 
(see Appendix A, part b). For 3 < 0, d^x is backward differenced to main- 
tain stability, and the first-order three-point backward difference is used. 

The scheme generalizes to three dimensions in the usual Douglas-Gunn-like 
fashion with the 6j^ operator replacing the usual identity operator. Finally, 
it is remarked that this procedure generalizes to any number of operators :in 
time that contain the factor (1 - E7^). For_^^xample, if, as in differencing 
(1), an operator such as (1 - E'^ ^).f (6x»6t)<t'j ,k appears, the operator £(6j^,6|.) 


replaces the usual identity operator of the Douglas-Gunn algorithm, 
stability is not generally ensured. 


Of course, 


c) A Numerical Verification of the Linear Algorithms 

The following example is presented as a simple check on the accuracy and 
stability of the implicit schemes. The motion is that of a parabolic arc 
airfoil increasing its thickness from zero to C.l, according to the relation 
shown in figure 2, in the time it takes a fluid particle to travel fifteen 
chord lengths at the free-stream velocity Since the motion is not 

periodic, t is scaled by c/U„ rather than 1/u). This is equivalent to 


A-6373 


7 


taking k “ I in (4) so that time is given in chord lengths traveled rather 
than radians. Solutions of the linear version of (4) (l.e. , y * -1) are com- 
pared in figure 2 in terms of pressure coefficient, Cp • at x ■ 0.525 

as computed from the second-order ADI scheme and from linear theory (exact) . 
Equivalent results were obtained with the second order AF scheme. Here the 
X grid spacing was uniform at Ax “ 0,05 chord lengths. The vertical mesh 
spacing, also uniform, was 0.1. The time step was five times greater than 
that permitted by the semi-implicit scheme (see (5)), and no accuracy or sta- 
bility problems were encountered. 


III. ONE-DIMENSIONAL SHOCK WAVE MOTION 


a) Model Problem 

The shock capturing properties of the implicit flnlte-diffeience schemes 
can be investigated using the invlscld Burger's equation, given here in terms 
of the velocity potential 




( 21 ) 


along with boundary and initial conditions 


ij)(0,t) = 0 , 


(Ji^(0,t) = 

L 


i}>(x,0) 


^ > 0 < X < Xg^ 

\ ■ W’‘ ■ "“o’ ■ 


Xg £ X i I 

o 


(22) 


The initial conditions are sketched In figure 3. The solution to (21) and (22) 

is (|ijj “ 4*x ^ ‘^’x ” “("x ^ ^ *s» where Xg is the 

L R 

instantaneous location of the shock wave and is given by 

Xg “ Xg + (here we exclude so-called "expansion shocks"). The 

potential i|i must be continuous for all x,t. 


b) Difference Schemes 

Equation (21) is written in difference form as 

* Vj - 0 C23) 

where 






) 


iua N 


'tuCOlILlTY OF THE 
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The quantity fj is some difference approximation for (<J)x)^i and the dif- 


ference operator remains to be specified. 


To avoid iteration at the n+1 
i.n+ ,.v,„ AP sclieme, 


level, fj is linearized by expandinj; in terras of . For the AF 

is evaluated at time level n+J , and the following first- and second-order 
linearizations can be made 


f p = = <)>" + 0(At) 

“l Xj ’Xj Xj 


(and by Taylor series expansion) 


AF 


° (•}>" + 241 " (4>”^* - 4>!! ) + O(At^) 


J 


j 




'j 


(24a) 


(24b) 


For the ADI scneme, 4>x averaged at the n and n+1 time levels, which, 
along with (24b) , gives 


^ADt 


j 


I [(♦ 


n+1 


) 




.n ,n+l 
p 4' 

"'j "'j 


0(At2) 


(24c) 


(The use of the terms "AF" and "ADI" is strictly not correct in this one- 
dimensional application, but the terminology is convenient in relatlon.stiip to 
the subsequent development of the two-dimensional algorithms.) 

In Section II it was determined that central and backward differences for 
4>jjjj were required to maintain stability for the cases (B > 0 (subsonic) and 
g < 0 (supersonic), respectively. In differencing (23) then, central and 
backward differences should be used for f when 4>x ^ (subsonic) and 

4>x > 0 (supersonic), respectively. As shown by Murman (ref. 9), care must be 
taken in switching from one difference operator to another, otherwise the con- 
servative form and its correct weak solution may not be maintained. Switched 
differences can be used while maintaining proper conservation form by express- 


ing 

D. 


in the form 


xfj = (<i - 






with £j = 0 
terminology 


j+(l/2) 

or 1 for (4>j+j “ less or greater than zero 


Ax 


(25) 


, n 


In Murman 's 


'j-1 

'j-1 


■j-1 


0 

1 

1 

0 


E . = 

3 

e . 


E . = 


E . = 


0 subsonic point 

1 => supersonic point 

0 sliock point 

1 =» parabolic (or sonic) point 


(26) 
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The procedure used here to maintain conservation form with switched dif- 
ferences can be summarized as follows: write the difference equation in con- 

servation form as in (23), linearize as <n (2A) , then apply the switching 
operator Dy^ as in (25). Violation of ,his procedure can result in noncon- 
aervative difference scliemes whicli yield erroneous shock speeds, as demon- 
strated in Appendix B. 


Using time linearizations consistent with the order of firot- and 

second-order accurate implicit differences for (21) are given (ay 


First order AF and ADI 
,n+l ,n+l 


♦p - ♦j':; “ 


3<P 


Second order ADI (two level) 
n+1 , , n+1 . , n+1 




Second order AF (three level) 




'j 




(27a) 


(27b) 


(27c) 


Note that the first-order AF and ADI schemes are equivalent. Also, the term 
"second order" Ln (27b) and (27c) applies only for locally subsonic regions 
because the upwind difference used in locally supersonic flows is 

only first-order accurate. 

For purposes of comparison, the semi- implicit scheme of reference 4 is 


,n+l 


- = -4xAtD^(<|)x^)2 + 


(28) 


and it is subject to the stability criterion At < 4x/(2| • 


c) Shock Characterization 

To investigate the shock-capturing characteristics of the numerical 
schemes, four types of shock motions are considered, as summarized in table I. 
The terms "subsonic" and "supersonic" are relative to the coordinate system 
and not to the moving shock. 
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TABLE I.“ CLASSIFICATION OF SHOCK MOTIONS 


Shock 

Cha'.acteristic 

) 

Spatial dif ferencingj 

1. 

Supersonic- to-supersonlc 

> 0, <!),, >0 
’'l ’‘r 

Backward 

1 

2. 

Subaonlc-to-subsonl( 

'^Xj^ '^Xj^ ® 

Central ! 

1 

3. 

Suporsonic-to-subsonic (down- 
stream moving) 

ij> +i|> <0 

^L % 

j 

Mixed 


Supersonic- to-subsonic 
(upstream moving) 

*L ^R 

Mixed 


A parameter that has a significant effect on the shock capturing proper- 
ties of any one of the schemes is the number of Ax increments the shock 
wave travels in a time At 


T 


At 

Ax dt 


(29) 


where the correct shock speed for (21) is 


dXg/dt » (1|)X^ + <t.Xj^), 


and a uniform 


gfld is assumed. A value ■> 1 corresponds to the case where the shock 
moves a distance of one v- jsSL^d point per time step. Another parameter of 
interest is the Courant number 


V 




(30) 


and the stability restriction for the semi- implicit scheme is |v| < 1 


d) Computed Results 

Now consider the implicit schemes (27) appH^d to the model problem (21) 
with initial conditions (22) for the first type of shock motion listed in 
table I, the supersonic-to-supersonlc case. In this case the term (4>x^)x 
approximated by backward differences throughout. Results of he first-order 
(27a), second-order AOI (27b), and second-order AF (27c) sche.nes arc compared 
in figure A. Results from each of the schemes are shown at two vaiue.s of T 
for the same time. The shock profiles shown remain essentially unchanged for 
all greater times. The exact solution is given by the dashed lines, and the 
shock location is indicated by Xg, 

Results for tl.t. second type of shock motion, the subsonic-to-subsonic 
case, are shown in figure 5. In this case the term ('S’x^)x approximatoii 
by central differences throughout. Again the shock profiles shown remain 
essentially the same for all larger n. The first-order results are similar 
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in appearance to those for zne supersonic-to-snpersor.lc case. The second- 
order results contain oscillations that increase with T. 


A sequence of shock profiles for an upstream-moving, supersonic-to- 
subsonic shock Is shown in figure 6. For the supersonic-to-subsonlc cases, 
mixed differences are used for (<|ix^)x* sequence of shock profiles is 

periodic, repeating every fifth time step, because T “ 0.2. All the schemes 
capture the shock sharply, but again oscillations are apparent in the second- 
order schemes. Similar results are shown in figure 7 for a downstream moving, 
supersonlc-to-subsonic shock. 


The first-order results of figure 7 are replotted in terms of i|) vs. x 
in figure 8, The dashed lines indicate the exact solution at different time 
steps. The exact shock location at time level n is indicated by the inter- 
section of two dashed lines and is marked by Xg on the abscissa. The loca- 

n 

tion of the shock relative to the mesh at n = 19 is repeated at n “ 24, and 
that at n “ 20 is repeated at n « 25, etc. Points identified as shock 
points ;'t level n according to (26) are denoted by Sj,. The solution to the 
difference equation (27a) for the case treated in rigure 3 is 


♦r‘ “ *i ' "“r. 


(31) 


where AXg 
At((j) 


the distance the shock travels in time At, is given by 
+ (j)jj ) . Hence, in updating from time level n to n+1, the 


AXg - 

solution remains unchanged for all (supersonic) poir' > to the left of the 


shock point. The shock point, and all (subsonic) points to the right of it, 
move to the dashed line that is the exact solution for n + 1. The test to 
determine assures that the shock point remains downstream of the shock 
for the value of T used here. 


e) Nonlinear Instab-'llty and Ovoi; shoots 

For large values of T, however, the shock-capturing procedure breaks 
down, and an instability occurs as illustrated in figure 9. Here the solution 
downstream of the exact shock location is correct for each n^ i , because of 
the test for the shock point can move only one grid point downstream per 

time step, while for T = 2, the ex, act shock location moves downstream at a 
rate of two grid points per time stt,.. An increasingly large discontinuity in 
develops, which appears ac a grov/ing overshoot in <|)x, and the process 
diverges. Reducing T to 1.1 slows the divergence rate, as illustrated in 
figure 10. For T = l-O. -be process is neutrally stable, as shown in 
figure 11, and the vs. x shock profile travels unattenuated after 

n = 3. For the case T "= 0*9, shown in figure 12, no instability occurs. The 
sequence is periodic with the solution repeating every tenth time step. An 
overshoot develops at n = 1 that decays with n until it repeats at n = 11. 


■;i .UlrlIUTY OF THE 
13 P(X)R 
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For this type of switched-differenclng-induced overshoot to occur, two 
conditions must be met: (1) the point immediately downstream of the shock 

must be a shock point, and (2) the (downstream moving) shock must move past 
this point in tim next time step. Referring to figure 3, these conditions can 
be expressed, respectively, as 


's+1 


or, equivalently. 


i|» < 0 X, - Xg > Ax 

^ B- 1 S ° 


■^X '^x , 

L R 


**'x ■ '^'x 

R 




X “X 

T(Ax/At) ^ 


‘*'x " K 


Ax 


< T 


(32a) 

(32b) 


(32c) 


This indicates that the occurrence of overshoots depends on the speed of the 
shock and its location relative to the mesh. As the shock propagates tlirough 
the grid, (32c) may be satisfied at some time levels and not at others. It 
follows from (32c) that these overshoots can never occur if 



This is a more restrictive condition than the Til requirement for stabil- 
ity. Inequalities can be derived from the upstream-moving shock case in a 
similar way. 

The type of instability and overshoot discussed here cannot occur in the 
seml-lraplicit scheme, because, for a downstream-moving supersonic- to-subsonic 
shock, the linear stability condif* n is more restrictive than the inequality 
(33). It is also interesting to note that for the nonconservative (in time) 
scheme detailed in the apper ,^x, no overshoots or instabilities occur. 

For the implicit schemes presented here, the shock-capturing process 
imposes time step limitations for both stability (due to the time lineariza- 
tion of the nonlinear term) and accuracy. The restriction T < 1 could be 
eliminated by fitting the shock as an internal boundary. Alternately, the 
followii simple correction process could perhaps be Implemented. First, 
update (p at n+1 by applying the implicit scheme as before. However, before 
proceeding to the next time step, any points crossed by the shock in time At 
would be corrected, i.e., they would be forced to lie on extrapolated curves 
from the upstream or downs -ream direction depending on the direction of shock 
travel. This procedure has been successfully applied in the one-dimensional 
case but has not yet been attempted in two-dimensions. 


IV. NUMERICAL ALGORITHMS FOR THE LOW-FREQUENCY TRANSONIC EQUATION 


The factored schemes introduced in Section II are combined here with the 
nonlinear mixed difference operators of Section III to form a complete algo- 
rithm for solving (4). 


a) Alternating Direction Algorithm 


An implicit differencing of the low-frequency, transonic, small- 
perturbation equation (4), using two time levels, can be denoted 


where 8 j 
6yy and 


6i<S 


Xt'' 


,n+l 


- 2 


+ 




(34) 


= 2kM„^, is defined by (40), and the difference operators for 
are specified by the relations 


6yy 


f - 1 

1 - 1 

' * 

y 

1 

2 

^k+l ^k 

^k - ViJ 

^k+1 " ^k-1 

V / 


(35a) 


"xt 


At 


7-1 


iVl " ""j-lj 


[2 + y(l - = (At)"l(l - E;1)6^ 


(35b) 


Here p = 0,1 for the first- and second-order differences, respectively. The 
unusual choice of (xj+i - Xj.^) used with the 6xt operator is required to 
maintain conservation form tor variable meshes. The operator (35b) is not the 
usual spatial second-order-accurate relation derived from a Taylor series 
expansion, but, for a sufficiently smooth grid variation, the difference is 
nevertheless second-order accurate. For example. Introduce C “ 4(x) with 
one to one mapping x = x(C) such that Cj - Cj-i = A? is constant over the 
grid, and for the three-point backward spatial operator 


(1 - e;;1)(3 - E~^) 

"'j + l ■ ""j-l 


<(> 


(1 - E-l)(3 - E-1)4. 

^ ^ + 0(A52)][e^ 


.1 1 f-'Tp 

J 2ACx^ + (A^)2 x^^^ + ... 


o(as2)] 

(36) 


where xj+i(0 is expanded about Xj(C) via Taylor series. 
In (34) 


where 


^ k " 1 - MS - 2 + l)Mco<}-x , Px 

L J.kJ j,l 


‘*’x. . “ ^‘*‘j + (l/2),k ' ‘^j-0/2),k^'^^^j + (l/2) " ^j-(l/2)^ 
J > ^ 


(37) 
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Linearizing as in Section III gives 


.n+1 

j.k 


+ 


f3f 1 

n 


j.k 


.n+1 

j.k 


- 


'i.kj 


( 38 ) 


Two-time level averaging then gives 

f 


j.k 


f 


, n+1 

J.k "'j^kj 


_ 1 v" 
- 2 


+ f" 1 

2l^j,k ^ 'j,kj 

+ (1 - m£)4>" 


j.k 


(39) 


where 


^j.k 


VV ,. = 1 - - (y + DkC'l' 


m.n 

X 


J.k 


The 


Djj operator is applied as In Section III 




+ e 


J-lp-(l/2),k ■ 'l-(3/2) 


j+(l/2),k 

.]) 


- f 


j-(l/2) 




(40) 


where ej = 0 or 1 for (Vj+(i/2),k "k ^)-(l/2),k) > 0 or < 0. (The same 
smooth grid restriction (36) applies to' the case ej =■ 1,) 

Incorporating the above differences in the ADI algorithm gives 

X sweep 

- ♦5_,) - 

y sweep 


(41) 


where fn k = .'1‘v- )• On the x sweep, a matrix is generated that is 

’ *j.k ^j,k 

lower tridiagonal for supersonic points (ej-i Cj = 1) and tridiagonal for 
subsonic points (^j_j = Ej = 0). For parabolic points (£j-i = 0, ej = 1) all 
the entries are zero on the row, and for shock points (tj_i = 1, tj =0), 
there are four nonzero entries on the row. For the present work, a quadra- 
diagonal solver was used that solves the equations like the Thomas algorithm 
for a tridiagonal matrix. On the y sweep, the only x differences are in 

which is backward differenced. Hence, the scheme is marched from upstream 
to downstream solving a tridiagonal matrix for each x = constant line of y 
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grid points. The test for ej is only vequired for the x sweep. Had the 
y sweep been taken first here, as in Section II, the test would have been 
required on both sweeps, thereby reducing the computer efficiency of the 
method. 


b) Approximate Factorization Algorithm 


For the AF scheme the conservativ time-space differencing was carried 
out in the transformed potential 


^ = - 62 X + (Y -t- 1)M^4. 

(42) 

so that the resulting equation 


3i^Xt = -(4*5/2)j^ + ^yy 

(43) 


has the same nonlinear term modeled in Section III. T.-je use of such a trans- 
form is unnecessary for the basic scheme described below but does lead to some 
computational simplifications. The transform also a'*jids a difficulty with 

* 2 

certain contrived time-linearizations such as ( 2 (t) 

' J > ^ ^ j 

in which the switching of difference operators ft r numerical =c4 
not correspond to the correct sonic line. For f ■'.ample, use of this lineariza- 
tion in (37) lends to a coefficient of chat changes sign at a differ- 
ent point than VP jik 

J « 

Variable grid differencing was implemented in the AF rcheme by coordinate 
stretching transforms, although the algorithm is given here in the untrans- 
forraed plane. Using the second-order time linearization (24b) and mixed dif- 
ferencing (25) described in Section III, the second-order AF algorithm is 
given by 

[I + AtF"^j^/(3M£)][6M2n-l(l - E-1) - 2AtC-^ 



= 2 M2n-l[(-5E;l + 6 - 


(2E-1 


- 3 + 




where 


+ 2AtF"^j^j[At|jl/(3M2)]6yy +| 




(44) 


t = t/k 




y 

X 



\ ■ 

‘\+i 

3-.’“ - 

h ■ 

‘S+i 


+ n ^)/2 
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“ "^j"l^j-l^^j“l.k ” ^j-2,k^*^x' 


and 


yy 




1 ) - 


n,^(l 


- E“l) 


The above algorithm only requires two double-indexed storage arr/iys with 
intermediate values held in ^^iTk* efficiency, the operator Fj should 

be simultaneously computed along rows for tlse rlght-liand side and the first 
inversion. In forming Fj along rows, the term ?j+i(^J+i,k ~ ^j,k) la com- 
puted and then shifted Lack to obtain ?j^'^j,k “ 'J'j-i k^ 

^j-l^^J~l*k “ ^j-2,k^* only a first-order-accurate differencing for 

is needed (i.e., (8a)), the simplified contrived time linearlxatior: (24n) 
should be used. 


V. SAMPLE iVO-DIMENStONAL COMPUTATION 

In this section results from the AUI and AF schemes are compared Cor the 
case of an airfoil executing unsteady motion characteristic of a helicopter 
rotor in forward flight. Particular attention is directed toward comparing 
the sJiock maturing properties of the first- and second-order methods Cor 
various tyj.es of shock motion in light of the one-dimensional results pre- 
sented in Section III, For the computations reported in this section, a uni- 
form X mesh spacing of Ax = 0.02 chord lengths was maintained over a dis- 
tance extending from one chord length upstream of the airfoil to t}jc airfoil 
trailing edge. Tlie x meslj spacing outside tills interval was stretched 
smoothly out to the boundaries, which were located more than thirty chord 
lengths from the airfoil. Tlie mesh spacing in the streara-normai direction was 
uniform at Ay = 0.04 (y unsealed), starting at the airfoil surface and 
extending for a distance of 0.2 chord _ 2 ngths. It was stretclied smoothly 
beyond this point out to the grid boundary, which was located a distance of 
fifty chord lengths from the airfoil. 

The motion of a blade element near Che tip of an advancing helicopter 
rotor can be simulated by an accelerating, niul, subsequently decelerating, 
airfoil (refs, 2 and 10). A motion that is simpler to treat computationally, 
but which produces a similar flowfleld, Is that of a thickening-thinning air- 
foil as sliown in figure 13. A parabolic-arc airfoil thickens from zero to 
0.1 in the time interval in which a fluid particle travels fifteen chord 
lengths (relative to the airfoil) at the free-stream velocity U„. The airfoil 
then thins to a flat plate after thirty cliord lengths of travel. (Here k = 1 
as in Section lie, so time is given in chord lengths of free-stream travel.) 
During the thickening pliase of Che motion, a shock wave forms and propagates 
downstream. As the airfoil thins, the shock reverses, moves upstream, and 
propagates off the front of the airfoil chasing the expansion \jave Chat 
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precedes it. The initial formation and downstrennt propagation of the shock 
wave is indicated by the pjessure distributions shown in figure 14. The 
second-order ADI and AF schemes show no plottable difference. 


An indication of the s!;ock speed can be obtained from the jump relations 
for (4). For a normal shock, the speed is given by 

(l/2k^^)[l - m£ + (1/2) (y + l)M!!!(Cp^ + (wliere Gp = -2<<.x) , and, in 


transonic small-disturbance theory, Cp* is given by setting V(. = 0 in (1) 


or (4), i.e., 1 - + (1/2) (y + DM^C 


* = 0. 


Hence, when the average 


across the shock is equal to Cp , the shock is stationary; otherwise the dif- 


ference between (C^ + C,j )/2 

R 


and 


indicates the direction and speed of 


the shock. Thus, at t => 11.5, the shock is traveling downstream, while at 
t “ 18.25, it is nearly stationary. As in Section III, this downstream-moving 
supersonic-to-subsonlc shock is properly captured, and the ADI and AF schemes 
give equivalent results. 


The subsequent upstream movement of the shock wave is shown in figure 15. 
At t - 26.875, there is oi't^y a small region of supersonic flow (for which 
Cp < Cp*), and the shock is traveling rapidly upstream. Here, as in Sec- 
tion III, the upstream-moving supersonlc-to-subsonlc shock \.rave produces over- 
shoots in the second-order scheme results. At later times, the shock is a 
subsonic-to-subsonic type, and it is accelerating in the upstream direction, 
as indicated by the increasingly large ratio of | (Cp + Cp )/2Cp*| us a 
function of Lime. Here, as in Section III, the f irst-order^^scheraes produce 
some "smearing" while the second-order schemes develop overshoots. Halving 
the time step reduced the amplitude of these oscillations in the second-order 
solut'ions and produced no plottable difference in the first-order solutions. 


Another view of the solution is afforded by the mid-chord pressure coeffi- 
cients as a function of time as shown in figure 16. There is a lag of about 
two chord lengths between the maximum flow expansion point and the maximum air- 
foil thickness. The effect of the shock wave propagating upstream past the 
mid-chord point is evident at about t = 26. 


It is interesting to compare the relative efficiencies of the implicit 
and semi-implicit schemes for the present example and to try to extrapolate the 
results of such a comparison to cases of more practical interest. In the 
present calculations, At was fixed at At = 0.125, where time is given here 
in chord lengths of free-stream travel relative to tlie airfoil. Considerably 
larger time steps could have been taken over much of the period of the motion, 
the major restriction being adequate resolution of the rapid shock motion 
upstream of Che airfoil. The time step restriction for Che sorai-impllcit 
scheme is given by 


At < Min . , , 

J 


1 - 


Ml + 


2m2Ax 

Y -f 1 


m"'c 

2 M^p 




(4 5) 


where the minimum is taken over all spatial grid points j,k in the flow- 
field. In the present calculations, Che most restrictive At according to 
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(45) would be At “ 0.0328 at t 16. Since the implicit schemes require 
about fifty percent more work per time step, the efficiency ratio of the 
implicit schemes relative to the semi- Implicit schemes 1s roughly 2.5 for this 
case. In more practical cases, i.e. , for blunt leading-edge airfoils at angle 
of attack (e.g. , see ref. 3), the leading edge singularity is much stronger, 
and a much finer Ax is required in the nose region. (The solution is not a 
valid approximation to the Euler equations in this region because of the vio- 
lation of the small-disturbance assumptions under which the theory is derived. 
However, experience Indicates that for sufficiently fine grids near the lead- 
ing edge, these inaccuracies are confined to a small region, whereas for 
coarser grids, they can affect the solution over the entire surface of the 
airfoil.) Small Ax and large Cp both contribute to a considerably more 
severe time-step restriction for such cases, significantly increasing the 
efficiency ratio in favor of the implicit schemes. 


CONCLUDING REMAEIKS 


Implicit two- and three-level approximate factorization finite-difference 
schemes for the low-frequency transonic small-disturbance equation have been 
constructed that have no time-step limitation based on a linear stability 
analysis. At this point neither scheme appea.- to be superior to the other 
either in terms of accuracy or efficiency, and the development and application 
of both will continue. Both schemes maintain stability for mixed flows by the 
use of conservatively-switched differences, i.e., central or backward differ- 
ences are used as the local flow condition is subsonic or supersonic. An 
instability, caused by the time linearization and switching test, occurs when- 
ever a shock wave travels more than one spatial grid point oer time step. A 
mild time-step restriction is thus imposed that is usually orders of magnitude 
less severe than the one associated with explicit schemes. The implicit 
schemes should prove useful in the analysis of unsteady airfoil motions, 
including pitch, plunge and free-stream Mach number oscillations; and tlie 
extension to three dimensions, as outlined in Section II, should present no 
additional difficulty. 
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APPENDIX A 


a) Stability of the AF Scheme 


To investigate the stability of uO) simplified boundary conditions are 
taken as either (a) spatial periodicity, or, because the all periodic case 
leads to a singular factor [(1 -• Ex^) “ bdyy], (b) spatial periodicity in x 
and specified function on the y boundaries. If ^ is the JxK vector 
ordered k ^ ° *'^*^*’ written 

= cr + ? 

The matrices A, B and C are either block diagonal or block tridiagonal with 
clrculant blocks (because of periodicity in x) and can be simultaneously 
diagonalized by unitary transforms. Hence, in an need only 

examine the eigenvtilues for stability. For the periodic boundary condition 
case, the differenc.- operators generate clrculant matrices with known eigen- 
values (refs. 11 and 12), and the matrix stability test is equivalent to the 
von Neumann test. Thus (10) is transformed 


(1 + a - ae^^J) 


in+l 

.k 


|l - e”^®j + 4C sin^l^jj^""^ 

= |l - e + ac(e ^ 


0. 

2| k 

2 


'4i 


j .k 


(Al) 


where c = b/(Ay)^ and 6j ZttJ/J, = 2Trk/K, or 0j^ - 2Trk/ (K + 1) if the 

function is specified at ^nd points jin y. This can be re^grouped as 


■I 


sin 


01 

2|_k 
2 


n+1 
k 


-iO 


Now 


= 1 - e 


;]j+ +«c 

+ ac^l - ♦" 


a^l - e - e = Aa sin^ 

so (A2) has the form 


2, 


land Ac sln2|_k 
2 ) 


> 1 


-n+l 


Since a > a > 0, |5j'|'kl > uncondlcional suability is assured, 

although in the periodic case when j = J and k = K, (A3) is formally 
indeterminate. 


a + Ic 
a + ic 


K a2 + 


imrn 




(A3) 


, r 0? 


OOB 
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b) Stability of the ADI scheme 


To demonstrate stability of ('8), assume a periodicity condition in x 
and Dirichlet conditions in y — however, assuming periodicity in x would 
make [6j^ - B(At/2)5j(j{] a singular operator , so its inversion in (18b) is not 
formally possible. Eliminating 'I’jTk (18b) using (18a) and the inversion 
of - (At/2)fiyy3 leads after some rearrangement to 


(6^ - == I'^x - T 




\ / 




(A4) 


Again the coefficient matrices can be simultaneously diagonalized by unitary 
transformations, and stability depends only on the eigenvalues. In the trans- 
formed variables, $ 



(A5) 


For sin^(0j/2) # 0 the left-hand coefficient can be divided out, and one can 
readil'. show that a ratio of eigenvalues appear which has the same type of 
inequa] cy as (A3). Thus for all roots such that sin^(0j/2) > 0, stability 
is unconditional. For the singular case, when sln(0j/2) = 0, the system is 
formally indeterminant. 
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APPENDIX U 


A scheme is considered here that Illustrates loss of conservation form by 
the use of proper conservative spatial differencing but with an improper time 
linearization. An interesting property of the scheme is that arbitrarily 
large time steps con be taken without generating t!ie type of instability dis- 
cussed in Section III. 


An implicit, first-order accurate approximation lo (21) is given by 


-1., . n+i n. 


which uses the notation defined in Section III and where . To avoid 

solving a nonlinear system at each time level, the first-order accurate 

/V n 

approximation (|ij - <Jij is made. Wltl) this approximation it can be shown by a 
flux summation argument that, for (Bl) applied to the model data shown li 
figure 3, the shock speed is 


dt 



(B2) 


whereas the correct shock speed for (21) is 


dt 



(B3) 


Thus (Bl) satisfies the correct shock conditions only for stationary shocks 

('t'x “ “‘f’x shocks with zero strength (ifij^ = <J>x ) » when At/Ax ->■ 0. 

L R L R 

in interesting feature of this nonconservative (in time) scheme is that 
as At/Ax », the shock speed approaches zero, as shown in figure 17. In 
fact, T -*■ (<l>x^ + 4>x 

bllities are expected to occur for large At, and numerical experiments con- 
firm this. This feature is attractive in applications in which only the 
steady-state solution is of interest, and using (Bl) with tlie first-order AF 
scheme does indeed provide a very effective relaxation algorithm. 




- ikv ) 

R 


as At/Ax ®. Hence, from (33), no insta- 
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igure la.- Characteristics for the sniall- l-igure lb*- Characteristics for the Xow- 

distivrbance equation (1). frequency approximation (4).. 
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(c) SECOND ORDER AF (THREE LEVEL) 


Figure 4.- Shock profiles for supersonic-to-supersonic case 
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(a) FIRST ORDER 


<t‘x 0 



(b) SECOND ORDER ADI (TWO LEVEL) 
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(c) SECOND ORDER AF (THREE LEVEL) 


Figure 7.- Shock profiles for downstream-moving, supersonic-to-subsonic case 





Figure 8.- vs x -ior downstream moving, supersonic-to-subsonic case, T = 0.2, v = 2.0 



Figure 9.- First-. Jer results for T = 2.0, v ~ 8.0. 










ADI SCHEME 

FIRST ORDER, At=0.125 Mqo = 0.85 

SECOND ORDER At ==0.125 m= 1.642 



Figure 15.- Pressure coefficient for upstream propagation of the shock wave. 



Figure 15.- Concluded. 
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Figure 16 ,- Midchord pressure coefficient vs time. 
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Figure 16.- Concluded. 
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